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Abstract 

Fundamental solitons pinned to the interface between three semi-infinite one-dimensional non- 
linear dynamical chains, coupled at a single site, are investigated. The light propagation in the 
respective system with the self-attractive on-site cubic nonlinearity, which can be implemented 
as an array of nonlinear optical waveguides, is modeled by the system of three discrete nonlin- 
ear Schrodinger equations. The formation, stability and dynamics of symmetric and asymmetric 
fundamental solitons centered at the interface are investigated analytically by means of the vari- 
ational approximation (VA) and in a numerical form. The VA predicts that two asymmetric and 
two antisymmetric branches exist in the entire parameter space, while four asymmetric modes 
and the symmetric one can be found below some critical value of the inter-lattice coupling param- 
eter - actually, past the symmetry-breaking bifurcation. At this bifurcation point, the symmetric 
branch is destabilized and two new asymmetric soliton branches appear, one stable and the other 
unstable. In this area, the antisymmetric branch changes its character, getting stabilized against 
oscillatory perturbations. In direct simulations, unstable symmetric modes radiate a part of their 
power, staying trapped around the interface. Highly unstable asymmetric modes transform into 
localized breathers traveling from the interface region across the lattice without significant power 
loss. 
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1. Introduction 

Surface modes, which are a special type of waves localized at interfaces between different 
media, were first predicted as localized Tamm electronic states at the edge of a truncated pe- 
riodic potential [1]. In optics, it was predicted theoretically and confirmed experimentally that 
the nonlinear self-trapping of light near the edge of a waveguide array with the self-focusing 
nonlinearity can lead to the formation of discrete surface solitons J2,|3t]. Various settings for 
the creation of surface solitons were also proposed for Bose-Einstein condensates (BECs) Jj]. 
The general framework for the description of such localized patterns is provided by the discrete 
nonlinear-Schrodinger (DNLS) equations, with appropriate boundary conditions [5D. 
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The solitary surface modes exist above a certain threshold power and, in a certain domain 
of the parameter space, different surface modes can exist simultaneously. The surface modes 
may be understood as discrete optical solitons localized near the surface |@]-|@], or as lattice 
solitons pinned by defects I^l- lfl4ll . Surface solitons supported by truncated superlattices were 
investigated too 111511 . Moreover, the light localization in self-defocusing nonlinear media in the 
form of surface gap solitons has been predicted and observed in Refs. HQ3, USB, and 



the concept of multi-gap surface solitons, i.e., mutually t rapped surface modes with components 
associated with different spectral gaps, was put forward | lj| 20tl_( multi-gap. alias "inter-gap", or 

"semi-gap", solitons are also known in uniform lattice media B21I1 ). A short review of surface 

i — i 

solitons in discrete systems was given in Ref. [22]. 

The studies of surface modes have shown that nonlinear discrete photonic and matter-wave 
systems support spatially localized states with sundry symmetries (which can be controlled by 
the insertion of suitable defects into the system) 11231 12411 . Related to this is the possibility of 
the spontaneous symmetry breaking (SSB) in symmetric dual-core systems, with a linear cou- 
pling between the two parallel cores. In fact, the SSB bifurcation, which destabilizes symmetric 
states and gives rise to asymmetric ones, was originally predicted in terms of the self-trapping 
in discrete systems B25I1 . In the physically important model of dual-core nonlinear optical fibers, 
the SSB instability was discovered in Ref. 12611 . and the respective bifurcations for continuous- 
wave states were studied in detail in Ref. [27], for various types of the intra-core nonlinearities. 
Further, the SSB was studied for solitons (rather than continuous waves) in the model of the 
dual-core fiber with the cubic (Kerr) nonlinearity |28, 2^1, for gap solitons in the models of dual- 
core ll30ll and tri-core lf3lll fiber Bragg gratings, and for matter-wave solitons in the BEC loaded 
into a dual-core potential trap, that may be combined with a longitudinal optical lattice 113211 . In 
addition, the SSB was also analyzed in models describing optical media with quadratic 13311 and 
cubic-quintic J34ll nonlinearities. 

The variational approximation (VA) IB5I1 makes it possible to study the SSB in dual-core 
systems in an analytical form 0281.13611 . It is relevant to mention that the VA may allow one not 
only to describe fundamental localized modes, but also correctly predict their stability [37]. 

Continuing the investigation of the surface fundamental modes in coupled one-dimensional 
(ID) lattice system |36[ 38|, in this work we study soliton complexes formed at the interface 
of three identical semi-infinite lattices which form a trilete configuration, see Fig. 1 below. 
Such lattices can be written in bulk silica by means of femtosecond optical pulses, which is the 
technique which has made the creation of other types of surface solitons possible [39]. A major 
motivation for studying such trilete systems is that they provide a fundamental setting for the 
study of the SSB which is different from the well-studied dual-core configurations 13 ll 14011 . 

The present system is modeled by three DNLS equations coupled at a single lattice site. In 
the framework of this system, we consider complexes built of three fundamental surface solitons, 
each carried by one of the three chains. Actually, these fundamental solitons are realizations of 
Tamm states in the present setting. We show that they form families of symmetric and asym- 
metric complexes, which exhibit the power-threshold behavior. Their stability and propagation 
dynamics are examined in detail, using the VA and numerical calculations. The existence regions 
for all families of the soliton complexes are produced by means of both approaches, the analyt- 
ical results being quite close to the numerical ones. Stability properties, predicted by dint of the 
numerically implemented linear-stability analysis, are verified by direct numerical simulations 
of the soliton evolution. 

The rest of the paper is structured as follows. The model is formulated in Section 2, which is 
followed by the consideration of the VA in Section 3. Numerical results are collected in Section 
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4, and the paper is concluded by Section 5. 



2. The model 

As said above, the trilete lattice configuration is formed by three identical semi-infinite sub- 
lattices linked at one site. The inter-site coupling constant is scaled to be C = 1 inside the lattices, 
while a different constant, e, accounts for the linkage between the lattices, having the same sign 
as C, see Fig. Q] Thus, the triple lattice is modeled by the following system of three coupled 
DNLS equations, 

d<pn 2 

i—r- + (<Pn+\ + <Pn-i) - (5n,o0n-i + sS„fi(tJ/ n + 6„) + \<p n \ (j)„ - 0, 
az 

i—r- + (iAh+i + ifrn-V - S n fi\j/ n -\ + s6 nfi ((f>„ + 0„) + \if/„\ if/„ = 0, (1) 
az 

d8 

l ~T + + Q n-\) - 5„fiG„-i + S5„fi{tp n + ij/n) + \6 n \ 2 0„ = 0, 

az 

where z is the propagation distance (assuming that the chains represent three semi-infinite arrays 
of optical waveguides), n — 0, 1 , is the discrete coordinate in the chain (N is the total number 
of sites in each lattice that was used in actual numerical calculations), S n $ is the Rronecker's 
symbol, and rescaling is used to make the on-site self-attraction coefficient equal to 1 . 




Figure 1 : Three coupled chains, linked by the modified linear coupling e, form the trilete lattice. 

Soliton solutions to Eqs. ([1) are looked for in the usual form, 

4> n = u„ exp (ifxz) , if/„ = v„ exp (ifxz) , 9„ - w„ exp (jyuz) , (2) 

where u„, v„ and w„ are real discrete functions, and jj is the propagation constant. The corre- 
sponding stationary equations, 



- fiu n + u n+ i + u n -\ - 6„flU n ^i + s{v„ + w n )5 nfi + u 3 n = 0, 
-/iv„ + v„+i + v„_i - 6 n flV„-i + s(u„ + w n )5„fi + v\ = 0, 

-yUW„ + W„+i + W„-\ - 5„fiW„-\ + s(u n + V n )6„fl + w\ = 0, 



can be derived from the Lagrangian, 



L = L u + L v + L w + 2s (uqvq + v wq + w u ) , 
3 



(3) 



(4) 



where the intrinsic Lagrangians of the three semi-infinite chains are, respectively, 



hi = V \-fwl + \ u t + 2u„u n+1 \, 

«=(> \ z / 



and the last term in Eq. © accounts for the coupling between them. 



3. The variational approximation 



The analytical approach to the study of the soliton solutions may be based on the VA (vari- 
ational approximation) 13511 . which was adapted to discrete systems in several earlier works 



1 36, 37, 38,|41|,|42|]. Following this method, we adopt the following ansatz: 

\u„, v„, w„) — {A, B, C)exp(-an), atn>0, (5) 
with amplitudes A, B and C treated as variational parameters. As concerns inverse width a, 



following Ref. [4311 we fix it through a solution of the linearized version of Eqs. (0 for "tails" of 



the discrete solitons at n — > oo, from where it follows 



s = e~ a = ft/2 - a/W2) 2 -1, (6) 

hence the solitons may exist for values of the propagation constant fi > 2. Note that Eq. © 
yields s < 1, which is essential for the validity of analytical results displayed below - see Eqs. 
© and (ED. 

The substitution of ansatz ([5]) into Eq. (0]i yields the corresponding effective Lagrangian, 

L = U +L2+L 3 + 2s(AB+AC + BC), (7) 

l-s 2 2(1 -s 4 ) 

t—u + 2s 1 A 

Li = B — — — + — -B 4 

l-s 2 2(1 -s 4 ) 

The Euler-Lagrange equations for amplitudes A, B and C are derived from this Lagrangian in the 
form of 

— f+2eB + C = 0, 
dA 

dL 2 

—^-+2s(A + C) = 0, 
dB 

^+2s(A + B) = 0, 
oC 
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or, in the explicit form, 

--A + — !— -A 3 +s(B + C) = 0, 
s 1 - s 4 

--B+ —^—B 3 +e(A + C) = 0, (8) 
s 1 - s 4 

--C + — !— -C 3 +s(A + B) = 0. 
s 1 - s 4 

These equations allow us to predict the existence of different symmetric and asymmetric 
interface modes. 

3.1. Existence regions for the interface solitons 

The solution for symmetric solitons, with A = B — C, is easily obtained fromEqs. ([8]): 

A = ± V(l - s 4 )^ 1 -2s). (9) 

The dependence of this amplitude on s for fixed fi = 5, i.e., s = ^5 — V2~T) /2 * 0.21, is plotted in 
Fig. Ufa). It follows from Eq. (0 and is shown in Fig. |2jb) that the existence of the VA-predicted 
symmetric solutions at a given value of /j is limited to 

e < s c = (2s)-' =L- ^ 2 -4) ■ (10) 




Figure 2: (color online) Results of the variational approximation for symmetric solitons. (a) The amplitude versus the 
inter-chain coupling constant, e, for a fixed propagation constant, /j = 5, i.e., s = 0.21 (the same value is used in all 
examples displayed below). In this panel, A2 = — A\, see Eq. (9). (b) In the (/.(, e) plane, the existence region for the 
fundamental symmetric solitons, marked by symbol "SyS", is e < £ c = l/(2.v). For /i = 5, e c « 2.4. 



The variational calculations have shown that the SyS complexes with very small power are 
created near the s c = l/(2s). In this parameter region, when the nonlinear interaction is negligi- 
ble, we can actually analyze the corresponding linear trilete lattice system. The straightforward 
analysis shows that only symmetric linear surface localized complexes with arbitrary amplitude 
can be formed in the linear trilete lattice system. These linear modes exist exactly at s — 1 /(2s), 
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which coincides with the critical value of the inter-lattice coupling constant (s c ) for existence 
of symmetric nonlinear modes. Therefore, we expect that the stable nonlinear SyS branch with 
small power emerges from the linear symmetric localized mode. 

The critical value of s at which the branch of asymmetric solitons bifurcates from the sym- 
metric family is predicted by solving Eqs. ([8]) for the soliton's amplitudes with infinitesimal 
differences between them: B = A + 6B, C — A + 5C. Straightforward algebraic manipulations, 
which include the linearization with respect to 5B and 6C, yield the value of s at the SSB bifur- 
cation point: Sb = 2/(7 s), see Fig. [3] Asymmetric solutions exist if the linear coupling is weaker 
than at the bifurcation point, i.e., at e < s\,. 




Figure 3: (color online) Black, red, and blue symbols display amplitudes A = C and B of different branches of asymmetric 
solutions of the isosceles type, which are predicted by the variational approximation. Only positive branches of B, 
but not their negative counterparts, are shown. The olive solid line denotes the antisymmetric ("AnS") branch, with 
A = 0, B = -C + 0, which corresponds to Eq. 1111 . while the dashed orange line denotes the amplitude of the symmetric 
solitons ("SyS")- Each solution has its negative counterpart (see the text). The critical values of the coupling constant, 
corresponding to this s, are displayed too by vertical lines, £b = 2/(7 s) and £ c = I /(2s). 

An analytical solution for asymmetric modes can be easily found in the particular case when 
one amplitude vanishes, A = 0. Then two other amplitudes are 

B - -C = + V(l -s 4 )0-' +e). (11) 

This solution may be naturally called antisymmetric. Note that, unlike symmetric solution ©, it 
exists for all values of s. In fact, the antisymmetric solution is identical to its counterpart found 
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in the two-chain system in Ref. rf42ll . although this does not mean that the stability properties of 
the antisymmetric solutions are identical in the two systems. 

Another analytical asymmetric solution (which may be called an isosceles mode) has A = 
C + B + 0. In this case, B can be eliminated, 

B ia = ~ [-A ± VA 2 - 4(A 2 - (s- 1 + e)(l - s 4 ))J , (12) 
while A has to be found as a numerical solution of the remaining equation, 

A3 + ( ! ~ (f ~ 1) A = *5 ( ! ~ '*) V-3A 2 +4(l-i 4 )(^ 1 (13) 

In accordance with the above result that Sb = 2 1 (Is) determines the location of the SSB bifurca- 
tion, six different branches of solutions to Eq. (TPTt (three with positive and three with negative 
B) exist in the area of s < 2/(7 s) of parameter space (s,fi), while only two branches (one with 
positive and one with negative B) are found at s > 2 /(7s), as shown in Fig. [3] for /i = 5. In 
this figure, only branches with A - C and B > are shown. It is worthy to mention that, in 
the region of < s < s\, in the figure, all the asymmetric branches, with numbers 1, 2, and 3, 
coexist with the symmetric branch and the antisymmetric one with A = 0; further, in the region 
of £b < e < £ c , only branch A\ coexists with the same pair, while in the region of s c < s only 
branches A 1 and A = exist. It is possible to prove that solutions of Eq. dS) with all the three 
amplitudes different, A + B + C, do not exist, which was also confirmed numerically. 

Finally, it is relevant to mention that the asymptotic form of the stationary solutions can be 
easily understood for both s — > and £ — > 00. Indeed, in the limit of s — > the system splits into 
three uncoupled semi-infinite lattices. In this case, for fixed propagation constant jj, one may 
have either the usual single-component surface soliton J2J], or the zero solution. Accordingly, 
Fig. [3] demonstrates that the amplitudes of all modes at e — take either a fixed value, which 
actually corresponding to the single semi-infinite lattice, or vanish. 

In the limit of s — > 00, a straightforward analysis of Eqs. © yields the following explicit 
solution, which represents the isosceles mode in this limit: 

A a W^fl-S 4 ),^ ±£^£(1 - s 4 ), 

where B * 1.138 is a root of equation {/3 2 f - 2 (p 2 f + 4B 2 - 4 = 0, and a = /3 3 /2 « 0.737. This 
expression is similar to the asymptotic limit of Eq. ( fTTT i. B = -C ~ ^b(1 - s 4 ). 

3.2. Stability of the fundamental solitons ( the Vakhitov-Kolokolov criterion) 

The stability of the discrete solitons predicted by the VA can be estimated by dint of the 
Vakhitov-Kolokolov (VK) criterion, according to which the necessary condition for the stability 
of fundamental solitons is dP/d/j > 0, where P is the total power (norm) of the soliton J4TL l44ll . 
The total power of the solutions corresponding to ansatz © is 

DO 

P = V (u\ + v 2 + w 2 ) = -(A 1 + B 2 + C\ (14) 

to l ~ s 

where the relation between s and fi is given by Eq. ((6). 
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For the symmetric solution (01 with A = B — C, dependence P(fj) is plotted in Fig. 3Ja) for 
fixed e = 0.5, 1 .5 and 3. The slope of the P(/S) curve is positive in the whole existence region of 
the symmetric solutions, which, according to the VK criterion, indicates the possibility for the 
existence of stable solutions. 

The P(ji) curve for the asymmetric solution branches is plotted in Fig. |4|b) for three different 
values of s. Only the asymmetric branch which corresponds to amplitudes {Ai,Bi} in Fig. [3] 
is VK-unstable in the narrow area of the existence region, at fi close to 2 in Fig. SJb). Other 
asymmetric solutions might be stable according to the VK criterion. 

For the antisymmetric solutions with A = and B - -C + 0, the power curves are presented 
in Fig. EJc). These solutions also might be stable according to the VK criterion in the whole 
existence region for s < 1. For higher values of e, a VK-unstable region is located near the 
above-mentioned existence border, /i = 2 [seeEq. ©]. 



□ e-].5 

<7 c=3 



(b) AS 




• £-=0.5 

□ E=1J 















Figure 4: The norm (power) P versus propagation constant /i for (a) symmetric solutions ("SyS"), (b) asymmetric 
(isosceles) solutions ("AS"), and (c) antisymmetric solutions, with A = 0, B = -C 4- ("AnS"). The corresponding 
values of the inter-lattice coupling constant, £, are indicated in the panels. Different curves labeled by indices 1, 2, 3 for 
fixed e in panel (b) correspond to solutions with the same numbers in Fig. [3] 



4. Numerical results 

The predictions of the VA were tested by solving stationary equations (0, using an algorithm 
based on the modified Powell minimization method 04511 . The initial guess for constructing 
fundamental solitons centered at the interface of three linked semi-infinite chains (Fig. [TJ was 
taken in the form of uq — vq — wq — A for symmetric modes or uq - A, vq - B, wq - C for 
asymmetric ones (with different A, B, C), while the amplitudes of the lattice field at other sites 
were set to be 0. Eventually, solutions different from symmetric ones were found solely with 
B = -C, A = 0, or with A — C + B, as predicted analytically. The results presented here are 
obtained for identical coupled chains of length N = 51. 

The stability of the stationary modes was first checked through the linear-stability analysis. 
As a result, the eigenvalue (EV) spectrum for modes of small perturbations was found, following 
the procedure developed in Refs. [38, 45]. The calculations were performed in parameter plane 
(e, fi). These results were verified by direct numerical simulations of the full system of equations 
([TJ- The simulations relied upon a numerical code which used the sixth-order Runge-Kutta algo- 



rithm, as in Refs. [38, 45]. The simulations were initialized by taking stationary soliton profiles, 
to which random perturbations were added. 

Typical shapes of symmetric and asymmetric solitons found in the numerical form are dis- 
played in Fig. [5] well fitting the corresponding VA predictions. The respective dependencies of 
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the solitons' amplitudes, A and B, on coupling constant e are displayed for all types of the soli- 
tons in Fig. 0(a), while the corresponding dependencies P{s) are shown in Fig. [7] From these 
figures it can be concluded that both the amplitude and power characteristics give qualitatively 
the same information about the localized surface modes. The numerical results demonstrate that 
the symmetric and three asymmetric branches coexist in certain bounded regions of the parameter 
space, similar to what was predicted in Fig. [3] 
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Figure 5: (Color online) Numerically obtained profiles of the fundamental interface solitons: (a) symmetric (e = 1.5), 
(b) asymmetric of type 1 (e = 1.5), (c) asymmetric of type 2 (e = 0.7), and (d) asymmetric of type 3 (e = 0.5) (these 
types are defined in Fig. [6}. Panel (a) depicts the numerically found SyS profile (solid line with symbols) along with 
its numerically predicted counterparts ("VA") (dashed line with symbols), while on panels (b)-(d) the numerically found 
solution profiles are shown by symbols and VA counterparts by solid lines. 

In general, the comparison of the numerical and variational results demonstrates that the pre- 
dictions of the VA for the existence region of the symmetric and asymmetric isosceles complexes 
of types 2 and 3 are very accurate. On the other hand, the numerically found existence regions 
for the antisymmetric solitons of type 3 and antisymmetric complexes are bounded, on the con- 
trary to the prediction of the VA. For example, the VA predicts that the antisymmetric complexes 
can exist for arbitrary e and fi > 2, while the numerical calculation shown the appearance of the 
upper limit with respect to e at fixed yu, see Fig. 0(b). The upper existence boundary for the 
antisymmetric solutions was found at extremely high values of the inter-chain coupling constant, 
s ~ 100. 
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Figure 6: (Color online) (a) The amplitudes of numerically generated symmetric solitons ("SyS", the orange dashed 
line), and asymmetric ("AS") ones, corresponding to branches (1), (2), and (3) - red (triangles), blue (circles) and black 
symbols (squares). These branches are counterparts of the branches with the same numbers, which were predicted by 
the variational approximation and displayed in Fig. [3] The solid olive line denotes the antisymmetric branch with 
A = 0, B = -C + 0, which also has its variational counterpart in Fig. [5] Another similarity to Fig. [J]is that each solution 
has its negative mirror image, only the solutions with positive amplitudes being displayed here. The dotted lines denote 
the critical values £b and e c obtained numerically, (b) The upper boundary of the existence region of the AnS and AS i 
complexes is plotted by the solid black line with circles and red dashed line with squares, respectively. The corresponding 
soliton complexes exist in the area below the boundary curves. 

The linear-stability analysis shows that a stability window exists for the symmetric soliton 
complexes in the region between e c [i.e., for jj close to fi c = l/(2s c ) + 2e c , as it follows from 
Eq. (TT~0b 1 and s\„ see Fig. HJa). These symmetric complexes are formed by solitons which are 
wider and possess smaller amplitudes than their counterparts belonging to unstable symmetric 
complexes. 

The dynamics of the stable symmetric soliton complex is illustrated by Fig. [9ja) which shows 
the evolution of its component. In the rest of their existence region, the symmetric complexes are 
unstable, their instability being accounted for by purely real EV pairs. Under small perturbations, 
an unstable symmetric soliton complex sheds off a part of its power and relaxes into a trapped 
interface asymmetric breathing complex with a smaller power. The dynamics of the component 
belonging to the unstable symmetric complex is shown in Fig. Hfb). 

The symmetric solitons with s = correspond to the usual surface solitons in uncoupled 
semi-infinite lattices. It is well known that such surface modes are stable in almost the whole 
existence region |2], which is provided by the balance between the interaction of the soliton with 
the surface and the bulk lattice. The instability of the symmetric soliton complex in the coupled 
lattice system is related to a stronger repulsion from the interface than the repulsion induced 
by the intra-lattice potential energy barrier far from critical e c , thus enforcing the perturbed 
strongly pinned soliton component solitons to shed off a part of their power. A consequence is 
the formation of more stable interface asymmetric breathing modes with a smaller power. 

The difference between the interface and intra-lattice potential energies may also explain an 
enhanced stability of symmetric complexes which are centered farther away from the interface. 
An example is plotted in Fig. [TO] where purely real EVs are presented for symmetric soliton 
complexes whose components are centered at distances n c = 1, 2, 3, 4 from the lattice interface 
in the corresponding lattices. Eventually, the symmetric complexes centered at n c > 4 are stable 
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Figure 7: (Color online) (a) The power of numerically generated symmetric solitons ("SyS", the orange dashed line), and 
asymmetric ("AS") ones, corresponding to branches (1), (2), and (3) - red (triangles), blue (circles) and black symbols 
(squares). The solid olive line denotes the antisymmetric branch with A = 0, B = -C ^ 0. The dotted lines denote the 
critical values £b and e c obtained numerically. 



in their entire existence region. We stress that only the soliton complexes of the symmetric type 
were found, trying to set the centers of the component solitons farther from the interface. 

The EV spectra for the asymmetric isosceles solitons, which are labeled by subscripts 1 and 
3 in Figs. [3]and|6] indicate the exponential instability in their entire existence regions, see Fig. 
Ob). This conclusion is confirmed by direct simulations, which show that perturbed asymmetric 
complexes of these types radiate a significant part of their power in the form of a breathing com- 
plex which moves across the corresponding chains, see Fig. |9| c )- The instability of asymmetric 
solutions can be associated with the trend of the system to relax toward an energetically prefer- 
able state in the presence of the two above-mentioned forces - the repulsion from the lattice 
interface, and the force induced by the bulk lattice, which is measured by the respective Peierls- 
Nabarro potential barrier J45ll . The antisymmetric solutions with A — 0, B — -C are shown to be 
unstable against oscillatory perturbations in the region of s < Sb, see Fig. |8jc), while a narrow 
region with exponentially unstable antisymmetric complexes is found near the upper boundary 
for their existence domain. Dynamical calculations confirm the predictions of the linear stability 
analysis. 

The isosceles soliton branch labeled by subscript 2 in Figs. [3] and [6] which is stable in the 
whole existence region according to the linear-stability analysis, is stable in direct simulations as 
well. Therefore, it can be concluded that the symmetry-breaking bifurcation at e = £b is related 
to the stability exchange between the destabilized symmetric complex and the two emerging 
isosceles asymmetric complexes, one of which is stable and the other exponentially unstable. In 
addition, the antisymmetric complexes change their stability in the neighborhood of the bifurca- 
tion point (from the oscillatory instability at s < s\, to the stability at e > Sb). 

Returning to the global existence diagrams, it is worthy to note that two areas with coexisting 
symmetric, asymmetric and antisymmetric solutions can be identified: the domain featuring the 
coexistence of stable symmetric, unstable asymmetric and stable antisymmetric complexes, and 
the domain where the symmetric, antisymmetric and two asymmetric complexes are unstable 
and one isosceles solution branch is stable. 
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Figure 8: (Color online) The purely real eigenvalues (EVs) versus the interface coupling constant e for fixed jti = 5 are 
shown in plots (a) and (b) for the symmetric and the asymmetric ("AS") complexes of types 1 and 3, respectively. The 
implication of the picture shown in panel (a) is that the symmetric solitons are stable in the interval of ej, < e < e c . In plot 
(c), the real part of the complex EV (red squares) and pure real EVs (black circles) are displayed for the antisymmetric 
solitons. Another asymmetric complex (the one with amplitudes At, in Figs.[5]and|6j is stable in the entire existence 
region, < £ < ej, . 




Figure 9: Typical examples of the evolution of components of perturbed three-soliton complexes: (a) a stable symmetric 
soliton with e = 1.69; (b) an unstable symmetric complex with e = 1.09; (c) a component of the unstable asymmetric 
complex of type 1 (see Fig. [3j with s = 1 .69. As well as in other figures, the propagation constant of unperturbed solitons 
here is fi = 5. 



As said above, the dynamics of the interface soliton complexes is strongly related to the 
balance of the interface and bulk-lattice potential energies. When the inter-lattice coupling is 
too small, the interface complexes of both the symmetric and asymmetric types formed by the 
fundamental solitons are unstable. The exception is (as actually mentioned above) the isosceles 
asymmetric complex with B > A — C (mode 2 in Fig. [3] and [6]), which is stable in its entire 
existence region. The strong repulsion from the interface, in comparison to the force from the 
bulk lattice, makes the formation of stable surface solitons difficult. By increasing the inter-lattice 
coupling, the energy difference between the interface and intra-lattice forces gets smaller, which 
makes it possible to create stable localized interface symmetric and antisymmetric complexes. 
With the further increase of the inter-lattice linkage against the intra-lattice coupling, which 
means proceeding to s > 2 for fixed ft, the highly unstable asymmetric complex and stable 
antisymmetric one are the only soliton species generated by the system. These findings are 
correlated with results reported for the two-chain version of the present system in Ref. l38ll . 
where the stable symmetric, antisymmetric, and antisymmetric branches have been found in a 
certain part of the corresponding existence regions. 
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Figure 10: (Color online) Purely real eigenvalues (EVs) versus the interface coupling constant e for fixed /i = 5. Funda- 
mental solitons are centered at n c = 1,2, 3,4 in each chain, forming the symmetric configuration. Dotted lines denote 
the boundary of the corresponding soliton-existence regions. It is obvious that the solitons shifted deeper into the lattices 
are more stable. 



5. Conclusion 

In this paper, we have analyzed the properties of fundamental localized interface complexes 
excited at the interconnection of three nonlinear semi-infinite chains. This system may be real- 
ized in nonlinear optics and BEC. In the framework of the system of three semi-infinite DNLS 
equations with the self-focusing on-site nonlinearity, coupled at the single site, we have found, 
by means of the VA (variational approximation) and numerical calculations, the threshold value 
of the inter-lattice coupling constant bounding the creation of symmetric surface soliton com- 
plexes. The existence of families of symmetric and asymmetric interface soliton complexes has 
been demonstrated. The variational predictions for the stationary modes were checked numeri- 
cally. The stability analysis shows that the symmetric soliton complexes, which are created as 
a stable solution branch at the critical value of the lattice inter-coupling parameter, destabilize 
at the bifurcation point, where two isosceles asymmetric soliton complexes are created, one of 
them stable and the other exponentially unstable. In addition, the stability of the antisymmetric 
complexes changes twice - at the bifurcation point and near the upper boundary of their exis- 
tence region. The third isosceles solution branch is exponentially unstable in its entire existence 
region. In other words, in the trilete system, the symmetric solution branch, three asymmetric 
branches and the antisymmetric one coexist in a part of the parameter space (past the symmetry- 
breaking bifurcation). These solution branches exist with their mirror-image counterparts. Direct 
simulations demonstrate that unstable symmetric complexes are transformed, radiating away a 
part of their power, into robust oscillating modes in the form of localized interface breathing 
complexes, while the unstable asymmetric complexes form breathers traveling away from the 



13 



interface through the lattice. The origin of this behavior of the surface modes is related to the 
balance between the repulsion from the surface and the intra-lattice interactions. 

Because the instability of asymmetric soliton complexes is stipulated by the repulsion from 
the interface (the linkage site, n = 0), it seems plausible that they may be stabilized by making 
the on-site nonlinearity stronger at this site. Another interesting possibility is to consider solitons 
of the vortex type (cf. vortices studied in other linearly-coupled tri-core systems ll3lll40ll '). The 
respective ansatz may be taken as \u„,v„,w„} = A {l, e m ^, e 2m / 3 J exp(-an), cf. ansatz Q. This 
vortex may be classified as one of the "off-site" type, as its virtual pivot is located between the 
sites. 
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